---
title: "Extended Cholesky decomposition in R"
author: "Stéphane Laurent"
date: "2017-12-03"
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    keep_md: no
tags: R, maths
highlighter: "pandoc-solarized"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE)
```

Let $S$ be a symmetric positive semidefinite matrix of order $d$ having rank $r$. 
An *extended Cholesky decomposition* of $S$ is a triplet $(L,M,P)$ consisting of a lower triangular $r\times r$-matrix $L$, a $(d-r) \times r$-matrix $M$, and a permutation matrix $P$ of order $d$ such that, setting 
$$
C = \begin{pmatrix} L & 0 \\ M & 0 \end{pmatrix},
$$
one has $PSP' = CC'$. Besides, setting
$$
\widetilde{C} = \begin{pmatrix} L & 0 \\ M & I_{d-r} \end{pmatrix},
$$
one has $S={(\widetilde{C}'P)}'I_d^r \widetilde{C}'P$ where $I_d^r$ is the $d\times d$-matrix $\begin{pmatrix} I_r & 0 \\ 0 & 0 \end{pmatrix}$.

The R function below calculates an extended Cholesky decomposition.

```{r, attr.source='.numberLines'}
extendedCholesky <- function(S){
  C <- suppressWarnings(chol(S, pivot=TRUE))
  d <- nrow(C)
  P <- matrix(0, d, d)
  P[cbind(1:d, attr(C,"pivot"))] <- 1
  r <- attr(C, "rank")
  return(list(L = t(C[seq_len(r), seq_len(r), drop=FALSE]), 
              M = t(C[seq_len(r), seq_len(d-r)+r, drop=FALSE]), 
              P = P))
}
```

Let's check:

```{r, attr.source='.numberLines'}
d <- 3
##~~ check for a rank 1 matrix ~~##
S <- tcrossprod(c(1:d))
#~ extended Cholesky of S ~#
EC <- extendedCholesky(S); P <- EC$P; L <- EC$L; M <- EC$M
#~ C matrix ~#
C <- cbind(rbind(L,M), matrix(0, d, d-ncol(L)))
all.equal(P %*% S %*% t(P), C%*%t(C))
#~ C tilde matrix ~#
Ctilde <- cbind(rbind(L,M), 
                rbind(matrix(0, nrow(L), d-nrow(L)), diag(d-nrow(L))))
all.equal(
  t(t(Ctilde)%*%P) %*% 
    diag(c(rep(1, nrow(L)), rep(0, d-nrow(L)))) %*% 
    (t(Ctilde)%*%P), 
  S)
##~~ check for a rank 2 matrix ~~##
S <- tcrossprod(c(1:d)) + tcrossprod(d:1)
#~ extended Cholesky of S ~#
EC <- extendedCholesky(S); P <- EC$P; L <- EC$L; M <- EC$M
#~ C matrix ~#
C <- cbind(rbind(L,M), matrix(0, d, d-ncol(L)))
all.equal(P %*% S %*% t(P), C%*%t(C))
#~ C tilde matrix ~#
Ctilde <- cbind(rbind(L,M), 
                rbind(matrix(0, nrow(L), d-nrow(L)), diag(d-nrow(L))))
all.equal(
  t(t(Ctilde)%*%P) %*% 
    diag(c(rep(1, nrow(L)), rep(0, d-nrow(L)))) %*% 
    (t(Ctilde)%*%P), 
  S)
##~~ check for a rank 3 matrix ~~##
S <- toeplitz(d:1)
#~ extended Cholesky of S ~#
EC <- extendedCholesky(S); P <- EC$P; L <- EC$L; M <- EC$M
#~ C matrix ~#
C <- cbind(rbind(L,M), matrix(0, d, d-ncol(L)))
all.equal(P %*% S %*% t(P), C%*%t(C))
#~ C tilde matrix ~#
Ctilde <- cbind(rbind(L,M), 
                rbind(matrix(0, nrow(L), d-nrow(L)), diag(d-nrow(L))))
all.equal(
  t(t(Ctilde)%*%P) %*% 
    diag(c(rep(1, nrow(L)), rep(0, d-nrow(L)))) %*% 
    (t(Ctilde)%*%P), 
  S)
```